Gas pipeline networks have tremendous economic importance. As of September 2016, there were more than 2,700,000 km of natural gas pipelines and more than 4,500 km of hydrogen pipelines worldwide. In the United States in 2015, natural gas delivered by pipeline networks accounted for 29% of total primary energy consumption in the country.
The invention relates to the control of a gas pipeline network for the production, transmission, and distribution of a gas. Examples of gas pipeline networks include 1) natural gas gathering, transmission, and distribution pipeline networks; 2) pipeline networks for the production, transmission, and distribution of hydrogen, carbon monoxide, or syngas; 3) pipeline networks for the production, transmission, and distribution of an atmospheric gas.
Gas pipeline networks have constraints on the rates of gas that can be produced, transmitted, and distributed. In some cases, the constraints associated with production of the gas ultimately limit the amount that can be transmitted and distributed. In other cases, constraints associated with compressing and transmitting the gas are more constraining than constraints on production rates. In all cases, it is desirable to maximize the capacity factor of a gas pipeline network (defined as the total quantity of gas supplied divided by the total capacity of the pipeline). Maximizing the capacity factor lowers the total unit cost associated with the production, transmission, and distribution of the gas, thereby ultimately lowering costs for both the operator of the gas pipeline network and the consumer of the gas.
Maximizing the capacity factor of a gas pipeline network is difficult for several reasons. First, there are constraints on pressures within the gas pipeline network, and the nonlinear relationship between flow and pressure drops makes it difficult to identify operating conditions that simultaneously satisfy pressure constraints while maximizing the capacity factor. The determination of these operating conditions that satisfy pressure constraints while maximizing the capacity factor must be done quickly and reliably.
Second, maximizing the capacity factor of a gas pipeline network is difficult because demand for a gas is dynamic and fluid. Customer consumption rates of the gas often change. Furthermore, the amount of gas that could potentially be received and used by the customer at any given time is not a single fixed quantity. This is especially true of customers for industrial gases such as oxygen and nitrogen. Often, an industrial gas customer is able to benefit from increased flows of a gas because these increased flows would enable the customer to increase the capacity or quality of their own manufacturing process. For example, petroleum refineries often receive hydrogen gas from a hydrogen gas pipeline network. It is often the case that receiving additional hydrogen gas would enable them to increase their production of desulfurized gasoline and diesel. A condition under which a customer would benefit from gas flows higher than their current consumption is referred to as latent demand for the gas.
The third difficulty in maximizing the capacity factor of a gas pipeline network is the information asymmetry between the operator of the gas pipeline network and the gas customers. Operators of the gas pipeline network are typically able to monitor conditions within the network to determine when there is additional capability to supply gas to a customer, whereas the customers typically do not have visibility to network operating conditions and constraints. On the other hand, while customers may know when they would benefit from increased flows of the gas, they may not always communicate this unmet demand to the operator of the gas pipeline network in an ongoing fashion. In short, customers typically lack information on unused gas supply capacity, and gas pipeline network operators typically lack information on latent customer demand. This information asymmetry tends to lower the capacity factor of the network. What is needed is a system and method for controlling a gas pipeline network to maximize the capacity factor, a system and method which simultaneously considers unused gas supply capacity and relevant latent customer demand.
A system and method for maximizing the capacity factor of a gas network provides set points to control elements which are operable to regulate pressure and flow. Control elements are operable to receive setpoints for the flow or pressure of gas at a certain location in the network, and use feedback control to approximately meet the setpoint. FIG. 1A illustrates an exemplary hydrogen gas pipeline network. This exemplary network illustrates at least certain of the physical elements that are controlled in accordance with embodiments of the present invention. Thus, control elements include pressure control elements 101 and flow control elements 102a, 102b. 
Industrial gas production plants associated with a gas pipeline network are control elements, because they are operable to regulate the pressure and flow of gas supplied into the network. Examples of industrial gas production plants include steam methane reformer plants 103 for the production of hydrogen, carbon monoxide, and/or syngas; and air separation units for the production of oxygen, nitrogen, and/or argon. These plants typically are equipped with a distributed control system and/or model predictive controller which is operable to regulate the flow of feedgas into the production plant and the flow and/or pressure of product gas supplied to the gas pipeline network.
Natural gas receipt points 104a, 104b are control elements, because they include a system of valves and/or compressors to regulate the flow of natural gas into the natural gas pipeline network.
Natural gas delivery points are control elements, because they include a system of valves and/or compressors to regulate the flow of natural gas out of the natural gas pipeline network.
Natural gas compressor stations are control elements, because they are operable to increase the pressure and regulate the flow of natural gas within a natural gas pipeline network.
Industrial gas customer receipt points 105 are control elements, because they are operable to receive a setpoint to regulate the flow and/or pressure of an industrial gas delivered to a customer.
Once available capacity and latent demand have been identified, setpoints can be received by flow control elements in order to increase the capacity factor of the network. To ensure that setpoints for flow control elements will result in satisfying demand and pressure constraints, it is necessary to calculate simultaneously the flows for each gas pipeline segment and gas pressures at network nodes. As described herein, in an exemplary embodiment, a network flow solution includes numerical values of flows for each pipeline segment and pressures for each pipeline junction that are: 1) self-consistent (in that laws of mass and momentum are satisfied), 2) satisfy customer demand constraints, and 3) satisfy pressure constraints.
The network flow solution may be determined using processing unit 110, an example of which is illustrated in FIG. 1B. Processing unit may be a server, or a series of servers, or form part of a server. Processing unit 110 comprises hardware, as described more fully herein, that is used in connection with executing software/computer programming code (i.e., computer readable instructions) to carry out the steps of the methods described herein. Processing unit 110 includes one or more processors 111. Processor 111 may be any type of processor, including but not limited to a special purpose or a general-purpose digital signal processor. Processor 111 may be connected to a communication infrastructure 116 (for example, a bus or network). Processing unit 110 also includes one or more memories 112, 113. Memory 112 may be random access memory (RAM). Memory 113 may include, for example, a hard disk drive and/or a removable storage drive, such as a floppy disk drive, a magnetic tape drive, or an optical disk drive, by way of example. Removable storage drive reads from and/or writes to a removable storage unit (e.g., a floppy disk, magnetic tape, optical disk, by way of example) as will be known to those skilled in the art. As will be understood by those skilled in the art, removable storage unit includes a computer usable storage medium having stored therein computer software and/or data. In alternative implementations, memory 113 may include other similar means for allowing computer programs or other instructions to be loaded into processing unit 110. Such means may include, for example, a removable storage unit and an interface. Examples of such means may include a removable memory chip (such as an EPROM, or PROM, or flash memory) and associated socket, and other removable storage units and interfaces which allow software and data to be transferred from removable storage unit to processing unit 110. Alternatively, the program may be executed and/or the data accessed from the removable storage unit, using the processor 111 of the processing unit 110. Computer system 111 may also include a communication interface 114. Communication interface 114 allows software and data to be transferred between processing unit 110 and external device(s) 115. Examples of communication interface 114 may include a modem, a network interface (such as an Ethernet card), and a communication port, by way of example. Software and data transferred via communication interface 114 are in the form of signals, which may be electronic, electromagnetic, optical, or other signals capable of being received by communication interface 114. These signals are provided to communication interface 114 via a communication path. Communication path carries signals and may be implemented using wire or cable, fiber optics, a phone line, a wireless link, a cellular phone link, a radio frequency link, or any other suitable communication channel, including a combination of the foregoing exemplary channels.
The terms “non-transitory computer readable medium”, “computer program medium” and “computer usable medium” are used generally to refer to media such as removable storage drive, a hard disk installed in hard disk drive, and non-transitory signals, as described herein. These computer program products are means for providing software to processing unit 110. However, these terms may also include signals (such as electrical, optical or electromagnetic signals) that embody the computer program disclosed herein. Computer programs are stored in memory 112 and/or memory 113. Computer programs may also be received via communication interface 114. Such computer programs, when executed, enable processing unit 110 to implement the present invention as discussed herein and may comprise, for example, model predictive controller software. Accordingly, such computer programs represent controllers of processing unit 110. Where the invention is implemented using software, the software may be stored in a computer program product and loaded into processing unit 110 sing removable storage drive, hard disk drive, or communication interface 114, to provide some examples.
External device(s) 115 may comprise one or more controllers that receive setpoint data from the software and are operable to control the network control elements described with reference to FIG. 1A.
It is difficult to calculate a network flow solution for a gas pipeline network because of a nonlinear equation that relates the decrease in pressure of a gas flowing through a pipeline segment (the “pressure drop”) to the flow rate of the gas.
This nonlinear relationship between flow and pressure drop requires that a nonconvex nonlinear optimization program be solved to calculate a network flow solution. Nonconvex nonlinear programs are known to be NP-complete (see Murty, K. G., & Kabadi, S. N. (1987). Some NP-complete problems in quadratic and nonlinear programming. Mathematical programming, 39(2), 117-129.). The time required to solve an NP-complete problem increases very quickly as the size of the problem grows. Currently, it is not known whether it is even possible to solve a large NP-complete quickly.
Embodiments of the present invention involve a system and method for controlling a gas pipeline network in order to maximize the capacity factor of the network. Embodiments of the invention determine whether it is hydraulically feasible to provide an increased flow rate of the gas to a customer. Embodiments of the invention further estimate whether the customer has latent demand for the gas. If it is hydraulically feasible to supply increased flow rates of the gas to the customer, and the customer has latent demand, then an updated request rate for the gas is received from the customer. In the preferred embodiment, the request from the customer is prompted by the gas production entity, who indicates to the customer that it is believed the customer has latent demand and the gas production entity has the capacity to supply the gas. The updated request rate is used to calculate a network flow solution, constituting flow rates for each pipeline segment and pressures for each pipeline junction. Elements of the network flow solution are received as setpoints by control elements.
Embodiments of the invention use a classification tree to determine whether a customer has latent demand for the gas.
Determining whether it is hydraulically feasible to supply an increased flow rate of the gas to a customer, as well as calculating a network flow solution, are enabled by several novel elements. First, the flow rate ranges for each pipeline segment are bounded under various scenarios for the supply and demand of the gas. These bounds are computed using a novel and computationally efficient network bisection method which is based on bounding the demand/supply imbalance on either side of a pipe segment of interest. Second, embodiments of the invention find the best linearization of the relationship between flow rate and pressure drop for each pipe segment, given the true nonlinear relationship between flow rate and pressure drop as well as the computed minimum and maximum flow rates for each segment. Third, embodiments of the invention use a linear program to compute a network flow solution, given the linearization of the relationship between flow rate and pressure drop for each segment. The linear program incorporates prior bounds on the inaccuracy of the pressure drop linearization to ensure that the network flow solution associated with an increased flow of gas to the customer will meet pressure constraints, given the actual nonlinear pressure drop relationship.
The following provides the notation used to describe the preferred embodiments of the invention. The first column identifies the mathematical notation, the second column is a description of the mathematical notation, and the third column indicates the units of measure that may be associated with the quantity.
Setsn ∈ NNodes (representing pipeline junctions)j ∈ AArcs (representing pipe segments and control elements)G = (N, A)Graph representing the layout of the gas pipeline networke ∈ {in, out}Arc endpoints(n, j) ∈ AinInlet of arc j intersects node n(n, j) ∈ AoutOutlet of arc j intersects node nn ∈ D ⊂ NDemand nodesn ∈ S ⊂ NSupply nodesj ∈ P ⊂ APipe arcsj ∈ C ⊂ AControl element arcsLj ∈ NLeft subgraph for arc jRj ∈ NRight subgraph for arc jParametersDjDiameter of pipe j[m]RGas constant[N m kmol−1 K−2]ZCompressibility factor[no units]LjLength of pipe j[m]MwMolecular weight of the gas[kg kmol−1]TrefReference temperature[K]ϵPipe roughness[m]αNonlinear pressure drop coefficient[Pa kg−1 m−1]fjFriction factor for pipe j[no units]μGas viscosity[Pa s]RejReynold's number for flow in pipe j[no units]qjminMinimum flow rate for flow in pipe j[kg/s]qjmaxMaximum flow rate for flow in pipe j[kg/s]bjIntercept for linear pressure drop model for pipe j[Pa2]mjSlope for linear pressure drop model for pipe j[Pa2 s/kg]ΔnMaximum additional amount to be supplied to[kg/s]customer in node nsnminMinimum production in node n[kg/s]snminMaximum production in node n[kg/s]VariablesdnDemand supplied rate in node n[kg/s]qjFlow rate in pipe j[kg/s]snProduction rate in node n[kg/s]pnnodePressure at node n[Pa]pjePressure at a particular end of a particular pipe[Pa]psnnodeSquared pressure at node n[Pa2]psjeSquared pressure at a particular end of a particular pipe[Pa2]psjerrMaximum absolute squared pressure drop error for pipe j[Pa2]psnerrMaximum absolute squared pressure error for node n[Pa2]
For the purposes of determining whether it is hydraulically feasible to provide an increased flow of gas to a customer, as well as for the purpose of computing a network flow solution, the layout of the pipeline network is represented by an undirected graph with a set of nodes (representing pipeline junctions) and arcs (representing pipeline segments and certain types of control elements). Here, some basic terminology associated with undirected graphs is introduced.
An undirected graph G=(N, A) is a set of nodes N and arcs A. The arc set A consists of unordered pairs of nodes. That is, an arc is a set {m, n}, where m,n ∈ N and m≠n. By convention, the notation (m, n), is used, rather than the notation {m, n}, and (m, n) and (n, m) are considered to be the same arc. If (m, n) is an arc in an undirected graph, it can be said that (m, n) is incident on nodes m and n. The degree of a node in an undirected graph is the number of arcs incident on it.
If (m, n) is an arc in a graph G=(N, A), it can be said that node m is adjacent to node n. The adjacency relation is symmetric for an undirected graph. If m is adjacent to n in a directed graph, sometimes it is written m→n.
A path of length k from a node m to a node m′ in a graph G=(N, A) is a sequence (n0, n1, n2, . . . , nk) of nodes such that m=n0, m′=nk, and (ni−1,ni) ∈ A for i=1,2, . . . , k. The length of the path is the number of arcs in the path. The path contains the nodes n0, n1, n2, . . . , nk and the arcs (n0, n1), (n1, n2), . . . , (nk−1, nk). (There is always a 0-length path from m to m). If there is path p from m to m′, it is said that m′ is reachable from m via p. A path is simple if all nodes in the path are distinct.
A subpath of path p=(n0, n1, n2, . . . , nk) is a contiguous subsequence of its nodes. That is, for any 0≤i≤j≤k, the subsequence of nodes (ni, ni+1, . . . , nj) is a subpath of p.
In an undirected graph, a path (n0, n1, n2, . . . , nk) forms a cycle if k≥3, n0=nk, and n1, n2, . . . , nk are distinct. A graph with no cycles is acyclic.
An undirected graph is connected if every pair of nodes is connected by a path. The connected components of a graph are the equivalence classes of nodes under the “is reachable from” relation. An undirected graph is connected if it has exactly one connected component, that is, if every node is reachable from every other node.
A graph G′=(N′,A′) is a subgraph of G=(N,A) if N′⊆N and A′⊆A. Given a set N′⊆N, the subgraph of G induced by N′ is the graph G′=(N′,A′), where A′={(m, n) ∈ A: m, n ∈ N′}.
To establish a sign convention for flow in a gas pipeline network represented by an undirected graph, it is necessary to designate one end of each pipe arc as an “inlet” and the other end as an “outlet”:
(n, j) ∈ AinInlet of arc j intersects node n(n, j) ∈ AoutOutlet of arc j intersects node n
This assignment can be done arbitrarily, as embodiments of the present invention allows for flow to travel in either direction. By convention, a flow has a positive sign if the gas is flowing from the “inlet” to the “outlet”, and the flow has a negative sign if the gas is flowing from the “outlet” to the “inlet”.
Some nodes in a network are associated with a supply for the gas and/or a demand for the gas. Nodes associated with the supply of a gas could correspond to steam methane reformers in a hydrogen network; air separation units in an atmospheric gas network; or gas wells or delivery points in a natural gas network. Nodes associated with a demand for the gas could correspond to refineries in a hydrogen network; factories in an atmospheric gas network; or receipt points in a natural gas network.
A set of mathematical equations govern flows and pressures within a gas pipeline network. These equations derive from basic physical principles of the conservation of mass and momentum. The mathematical constraints associated with a network flow solution are described below.
Node mass balance
The node mass balance stipulates that the total mass flow leaving a particular node is equal to the total mass flow entering that node.
            d      n        +                  ∑                  j          |                                    (                              n                ,                j                            )                        ∈                          A              in                                          ⁢                          ⁢              q        j              =                    ∑                  j          |                                    (                              n                ,                j                            )                        ∈                          A              out                                          ⁢                          ⁢              q        j              +          s      n      
The left-hand side of the equation represents the flow leaving a node, as dn is the customer demand associated with the node. The term Σj|(n,) ∈ Ain qj represents the flow associated with pipes whose “inlet” side is connected to the node. If the flow qj is positive, then it represents a flow leaving the node. The right-hand side of the equation represents the flow entering a node, as sn is the plant supply associated with the node. The term Σj|(n,j) ∈ Aout qj represents the flow associated with pipe segments whose “outlet” side is connected to the node. If the flow term qj is positive, then it represents a flow entering the node.
Node pressure continuity
The node pressure continuity equations require that the pressure at the pipe ends which is connected to a node should be the same as the pressure of the node.pjin=pnnode∀(n,j)∈Ain pjout=pnnode∀(n,j)∈Aout 
Pipe pressure drop
The relationship between the flow of a gas in the pipe is nonlinear. A commonly used equation representing the nonlinear pressure drop relationship for gas pipelines is presented here. Other nonlinear relationships may be used in connection with alternative embodiments of the present invention.
This nonlinear pressure drop equation for gases in cylindrical pipelines is derived based on two assumptions. First, it is assumed that the gas in the pipeline network is isothermal (the same temperature throughout). This is a reasonable assumption because pipelines are often buried underground and there is excellent heat transfer between the pipeline and the ground. Under the isothermal assumption, an energy balance on the gas in the pipeline yields the following eauation:
                    (                  p          j          in                )            2        -                  (                  p          j          out                )            2        =            q      j        ⁢                        q        j                    ⁢                            4          ⁢                                          ⁢          ZRT                                      M            w                    ⁢                      π            2                    ⁢                      D            j            4                              ⁡              [                                            4              ⁢                                                          ⁢                              f                j                            ⁢                              L                j                                                    D              j                                +                      2            ⁢                                                  ⁢                          ln              ⁡                              (                                                      p                    j                    in                                                        p                    j                    out                                                  )                                                    ]            
For gas pipelines, because the pipe lengths are large relative to the diameters, the term
      4    ⁢                  ⁢          f      j        ⁢          L      j        Dis so much greater than the term
  2  ⁢          ⁢      ln    ⁡          (                        p          j          in                          p          j          out                    )      that the latter term can be neglected. Under this assumption, then the nonlinear pressure drop relationship reduces to:(pjin)2−(pjout)2=αqj|qj|
with
  α  =            16      ⁢                          ⁢              ZRf        j            ⁢              T        ref            ⁢              L        j                            M        w            ⁢              π        2            ⁢              D        j        5            
where Z is the compressibility factor for the gas, which in most pipelines can be assumed to be a constant near 1; R is the universal gas constant; Tref is the reference temperature; Lj is the length of the pipeline segment; and the term fje is a friction factor for a pipe segment, which varies weakly based on the Reynolds number of flow in the pipe, and for most gas pipelines is in the range 0.01-0.08. Below is provided an explicit formula for the friction factor in terms of the Reynold's number. The dimensionless Reynold's number is defined as
            Re      j        =                  4        ⁢                                        q            j                                                π        ⁢                                  ⁢                  D          j                ⁢        μ              ,where μ is the gas viscosity.
If the flow is laminar (Reje<2100) then the friction factor is
      f          j      ,      L        =      64          Re      j      
If the flow is turbulent (Reje>4000), then the friction factor may be determined using the implicit Colebrook and White equation:
      1                  f                  j          ,          TR                      =            -      2        ⁢                  ⁢                            log          10                ⁡                  (                                    ϵ                              3.71                ⁢                                                                  ⁢                D                                      +                          2.51                                                Re                  j                                ⁢                                                      f                    j                                                                                )                    .      
An explicit expression for the friction factor for turbulent flow that is equivalent to the
Colebrook and White equation is
            f              j        ,        TR              =          1                        [                                    c              ⁡                              [                                                      W                    0                                    ⁡                                      (                                                                  e                                                  a                          /                          bc                                                                    /                      bc                                        )                                                  ]                                      -                          a              /              b                                ]                2              where            a      =              ϵ                  3.71          ⁢                                          ⁢          D                      ,          b      =              2.51        Re              ,                  and        ⁢                                  ⁢        c            =                        2                      ln            ⁡                          (              10              )                                      =        0.868589            and W0(·) is the principal Lambert-W function. See (More, A. A. (2006). Analytical solutions for the Colebrook and White equation and for pressure drop in ideal gas flow in pipes. Chemical engineering science, 61(16), 5515-5519) and (Brkic, D. (2009). Lambert W-function in hydraulics problems. In MASSEE International Congress on Mathematics MICOM, Ohrid.)
When the Reynolds number is between 2100 and 4000, the flow is in a transition range between laminar and turbulent flow and the accepted approach in the literature is to interpolate the friction factor between the laminar and the turbulent value, based on the Reynolds number, as follows:fj,TS=fj,\2100β+fj,TF|4000(1−β)with β=(4000-Rej)/(4000−2100).
Calculating Whether it is Hydraulically Feasible to Supply Additional Product to A Customer
A key enabler for controlling a gas pipeline network to increase its capacity factor is to determine whether it is hydraulically feasible to supply an increased flow rate of the gas to customer n, the increased flow rate being as much as Δn greater than the current flow rate of gas to the customer.
Because the relationship between pressure drop and flow is highly nonlinear, and because it is an NP-complete problem to determine the feasibility of supplying additional product in a gas pipeline network using this nonlinear pressure drop relationship, described is a method for determining the hydraulic feasibility of supplying additional product using a linearized pressure drop model.
FIG. 2 illustrates the nonlinear relationship between pressure drop and flow. The true nonlinear relationship is indicated by the solid line. If one approximates the true nonlinear relationship with a linear fit centered around zero, the linear fit severely underestimates the pressure drop for flow magnitudes exceeding 20. If one does a linear fit of the true pressure drop relationship in the range of flows between 15 and 20, the quality of the pressure drop estimate for negative flows is very poor. If one does a linear fit of the true pressure drop relationship in the range between −20 and −15 MMSCFD, the pressure drop estimate for positive flows is very poor.
To produce an accurate linearization of the pressure drop relationship for pipe segments, it is critical to bound the range of flow rates for each pipe segment. In examples below, a linearization based on tightly bounded flow rates is called a “tight linearization”. But note that if a customer receives additional product, this could alter the flow rates in the pipeline network.
Bounds on flow rates for a range of flow scenarios can be determined using mass balances and bounds on production for plants and demands from customers, even in the absence of any assumptions about pressure constraints and pressure drop relationships.
One method for bounding flows in pipeline segments based on mass balances is to formulate and solve a number of linear programs. For each pipe segment, one linear program can be used to determine the minimum flow rate in that segment and another linear program can be used to determine the maximum flow rate in that segment.
An exemplary embodiment of the present invention involves a method of bounding the flow rate in pipeline segments, under a range of demand/supply scenarios including the scenario where customer n takes additional product in a quantity up to Δn. The novel method is simple and computationally more efficient than the linear programming method.
For the pipe segment of interest (assumed to not be in a graph cycle), the pipeline network is bisected into two subgraphs at the pipe segment of interest: a “left” subgraph and a “right” subgraph associated with that pipe. Formally, the left subgraph Lj associated with pipe j is the set of nodes and arcs that are connected with the inlet node of pipe j once the arc representing pipe j is removed from the network. Formally, the right subgraph Rj associated with pipe j is the set of nodes and arcs that are connected with the outlet node of pipe j once the arc representing pipe j is removed from the network. Given the bisection of the flow network into a left subgraph and a right subgraph, it is then possible to calculate the minimum and maximum signed flow through pipe segment j, based on potential extremes in supply and demand imbalance in the left subgraph and the right subgraph.
To bound the flow rate in each pipeline segment, some quantities describing the imbalance between supply and demand are defined in the left and right subgraphs. The minimum undersupply in the left subgraph for pipe j is defined as sLjmin=(Σn ∈Lsnmin)−(Σn ∈Ldn+Δn). The minimum unmet demand in the right subgraph for pipe j is defined as dRjmin=(Σn ∈Rsnmax). The maximum oversupply in the left subgraph for pipe j is defined as sLmax=(Σn ∈Lsnmax)−(Σn ∈L dn). The maximum unmet demand in the right subgraph for pipe j is defined as dRmax=(Σn ∈Rdn+Δnmin).
Given the definitions above, the minimum and maximum feasible signed flow in the pipe segment are given by:qjmin=max {sLimin,dRimin },qjmax=min {sLimax,dRimax},
The equation for qjmin indicates that this minimum (or most negative) rate is the maximum of the minimum undersupply in the left subgraph and the minimum unmet demand in the right subgraph. The equation for qjmax indicates that this maximum (or most positive) rate is the minimum of the maximum oversupply in the left subgraph and the maximum unmet demand in the right subgraph.
The equations in the previous paragraph for calculating qjmin and qjmax can be derived from the node mass balance relationship, as follows. The node mass balance relationship, which was previously introduced, is
            d      n        +                  ∑                  j          |                                    (                              n                ,                j                            )                        ∈                          A              in                                          ⁢                          ⁢              q        j              =                    ∑                  j          |                                    (                              n                ,                j                            )                        ∈                          A              out                                          ⁢                          ⁢              q        j              +                  s        n            .      
Consider the left subgraph associated with pipe j. The left subgraph contains the node connected to the inlet of pipe j. Consider collapsing the entire left subgraph into the single node connected to the inlet of pipe j. Then,
      q    j    in    =                    ∑                  n          ∈                      L            j                              ⁢                          ⁢              s        n              -          d      n      
An upper bound for the inlet flow is qjin≤Σn ∈Ljsnmax−dn, and a lower bound for the inlet flow is qjin≥Σn ∈Ljsnmin−(dn+Δn). Similarly, an upper bound for the outlet flow is qjout≤Σn ∈R(dn+αn)−snmin and a lower bound is qjout≤Σn ∈Rjdn−snmax.
At steady state, the pipe inlet flow equals the outlet flow and
                              ∑                      n            ∈                          L              j                                      ⁢                                  ⁢                  s          n          min                    -              (                              d            n                    +                      Δ            n                          )              ≤                            ∑                      n            ∈                          R              j                                      ⁢                                  ⁢                  d          n                    -              s        n        max              ≤          q      j      in        =            q      j      out        =                  q        j            ≤                                    ∑                          n              ∈                              R                j                                              ⁢                                          ⁢                      (                                          d                n                            +                              Δ                n                                      )                          -                  s          n          min                    ≤                                    ∑                          n              ∈                              L                j                                              ⁢                                          ⁢                      s            n            max                          -                              d            n                    .                    Equivalently,
                    max        ⁢                  {                                                                      ∑                                      n                    ∈                                          L                      j                                                                      ⁢                                                                  ⁢                                  s                  n                  min                                            -                              (                                                      d                    n                                    +                                      Δ                    n                                                  )                                      ,                                                            ∑                                      n                    ∈                                          R                      j                                                                      ⁢                                                                  ⁢                                  d                  n                                            -                              s                n                max                                              }                    ≤              q        j        in              =                  q        j        out            =                        q          j                ≤                  min          ⁢                      {                                                                                ∑                                          n                      ∈                                              R                        j                                                                              ⁢                                                                          ⁢                                      (                                                                  d                        n                                            +                                              Δ                        n                                                              )                                                  -                                  s                  n                  min                                            ,                                                                    ∑                                          n                      ∈                                              L                        j                                                                              ⁢                                                                          ⁢                                      s                    n                    max                                                  -                                  d                  n                                                      }                                                  ⁢    or                      ⁢                            q          j          min                =                                            max              ⁢                              {                                                      s                                          L                      j                                        min                                    ,                                      d                                          R                      j                                        min                                                  }                                      ≤                    =                                                    q                j                            ≤                              min                ⁢                                  {                                                            s                                              L                        j                                            max                                        ,                                          d                                              R                        j                                            max                                                        }                                                      =                          q              j              max                                          ,      which completes the proof.
The bisection method for bounding flow rates in pipe segments is illustrated with an example. An example flow network is depicted in FIG. 3. This flow network has four customer demand nodes (nodes 1, 9, 12, and 16), and four plant supply nodes (nodes 2, 10, 13, and 17). This particular example relates to determining whether it is hydraulically feasible to supply additional product to the customer located at node 1. In this case, the current flow rate of the industrial gas to the customer at node 1 is 9 kg/s. It is determined whether it is possible to supply up to an additional 10 kg/s of gas flow rate, for a new total of 19 kg/s.
In FIG. 4, bounding the flow rate in the pipe going from node 5 to the customer at node 1 is illustrated. Per the indicated sign convention, a flow rate is negative if flow is in the direction leading from a higher-numbered node to a lower number node. The results of the bisection method indicate that the minimum signed flow rate in this pipe is −17 kg/s (corresponding to a flow of 17 kg/s to the customer at node 1), whereas the maximum signed flow rate in this pipe is −9 kg/s (corresponding to a flow of 9 kg/s to the customer at node 1). Thus, given plant production constraints and demand elsewhere in the network, we can only supply between 9 and 17 kg/s to the customer at node 1. This is less than the 19 kg/s originally envisioned, but is still significantly more than the current rate of 9 kg/s.
In FIG. 5, bounding the flow rate in the pipe going from the plant at node 10 to the junction at node 11 is illustrated. The results indicate that the minimum and maximum flows in this pipe are 7 kg/s and 12 kg/s, respectively, which is exactly consistent with the minimum and maximum production rate of the plant at node 10.
In FIG. 6, bounding the flow rate in the pipe going from node 3 to node 15 is illustrated. The results indicate that the minimum and maximum signed flows in this pipe are −6 kg/s and 5 kg/s, respectively. This indicates that flow can potentially go in either direction in this pipe.
FIG. 7, which shows data from computational experiments performed using Matlab on a computer with an Intel Core I 2.80 GHz processor, shows that the network bisection method for bounding the flow in pipeline segments is between 10 and 100 times faster than the linear programming method.
Finding the Best Linear Pressure-Drop Model Given a Scenario of Increased Customer Flow
The next step in assessing the hydraulic feasibility of providing additional flow to a customer is to linearize the nonlinear pressure drop relationship for each pipe, based on the flow bounds established for each pipe. This can be done analytically (if the bounded flow range is narrow enough that the friction factor can be assumed to be constant over the flow range), or numerically (if the bounded flow range is sufficiently wide that the friction factor varies significantly over the flow range). The sections below describe how a linearization can be accomplished either analytically or numerically. A linear pressure drop model is sought of the formpsjin−psjout=mjqj+bj∀j ∈P. 
Note that the fact that the flow range is bounded is critical to produce a good linear model. Without these bounds, a naïve linear model may be produced which is based on linearizing the nonlinear relationship about zero with a minimum and maximum flow magnitude equal to the total network demand. As will be shown in examples below, this generally does not produce good network flow solutions.
Finding the Least-Squares Linear Pressure-Drop Model Analytically: Slope-Intercept Form
If the bounded flow range is fairly narrow, then the friction factor as well as the nonlinear pressure drop coefficient a will be nearly constant and an analytical solution may be found for the least squares linear fit of the nonlinear pressure drop relationship.
Least squares solution for a linear model with g=qrmin and h=qjmax 
      (                  m        j        *            ,              b        j        *              )    =      arg    ⁢                  min                  m          ,          b                    ⁢                        ∫          g          h                ⁢                                            (                                                α                  ⁢                                                                          ⁢                  q                  ⁢                                                          q                                                                      -                mq                -                b                            )                        2                    ⁢          dq                    Evaluating the definite integral:
            ∫      g      h        ⁢                            (                                    α              ⁢                                                          ⁢              q              ⁢                                              q                                                      -            mq            -            b                    )                2            ⁢      dq        =                    b        2            ⁢      h        -                  b        2            ⁢      g        -                  g        3            ⁡              (                                            m              2                        3                    -                                    2              ⁢              α              ⁢                                                          ⁢              b              ⁢                                                          ⁢                              sign                ⁡                                  (                  g                  )                                                      3                          )              +                  h        3            ⁡              (                                            m              2                        3                    -                                    2              ⁢              α              ⁢                                                          ⁢              b              ⁢                                                          ⁢                              sign                ⁡                                  (                  h                  )                                                      3                          )              -                            α          2                ⁢                  g          5                ⁢                              sign            ⁡                          (              g              )                                2                    5        +                            α          2                ⁢                  h          5                ⁢                              sign            ⁡                          (              h              )                                2                    5        -                  bg        2            ⁢      m        +                  bh        2            ⁢      m        +                  α        ⁢                                  ⁢                  g          4                ⁢        m        ⁢                                  ⁢                  sign          ⁡                      (            g            )                              2        -                  α        ⁢                                  ⁢                  h          4                ⁢        m        ⁢                                  ⁢                  sign          ⁡                      (            h            )                              2      
This quantity is minimized when the partial derivatives with respect to b and m are simultaneously zero. These partial derivatives are
                    ∂                              ∫            g            h                    ⁢                                                    (                                                      α                    ⁢                                                                                  ⁢                    q                    ⁢                                                                q                                                                              -                  mq                  -                  b                                )                            2                        ⁢            dq                                      ∂        b              =                  2        ⁢                                  ⁢        bh            -              2        ⁢                                  ⁢        bg            -                        g          2                ⁢        m            +                        h          2                ⁢        m            +                        2          ⁢                                          ⁢                      ag            3                    ⁢                      sign            ⁡                          (              g              )                                      3            -                        2          ⁢                                          ⁢                      ah            3                    ⁢                      sign            ⁡                          (              h              )                                      3                                ∂                              ∫            g            h                    ⁢                                                    (                                                      α                    ⁢                                                                                  ⁢                    q                    ⁢                                                                q                                                                              -                  mq                  -                  b                                )                            2                        ⁢            dq                                      ∂        m              =                  bh        2            -              bg        2            -                        2          ⁢                                          ⁢                      g            3                    ⁢          m                3            +                        2          ⁢                                          ⁢                      h            3                    ⁢          m                3            +                                    ag            4                    ⁢                      sign            ⁡                          (              g              )                                      2            -                                    ah            4                    ⁢                      sign            ⁡                          (              h              )                                      2            
Setting the partial derivatives equal to zero, and solving for b and m, it is found that the form of the slope-intercept least squares linear model is:
            b      *        =          -                                                                  (                                                      α                    ⁢                                                                                  ⁢                                          g                      5                                        ⁢                                          sign                      ⁡                                              (                        g                        )                                                                              -                                      α                    ⁢                                                                                  ⁢                                          h                      5                                        ⁢                                          sign                      ⁡                                              (                        h                        )                                                                              -                                      8                    ⁢                                                                                  ⁢                    α                    ⁢                                                                                  ⁢                                          g                      3                                        ⁢                                          h                      2                                        ⁢                                          sign                      ⁡                                              (                        g                        )                                                                              +                                      8                    ⁢                                                                                  ⁢                    α                    ⁢                                                                                  ⁢                                          g                      2                                        ⁢                                          h                      3                                        ⁢                                          sign                      ⁡                                              (                        h                        )                                                                              +                                                                                                                                              α                    ⁢                                                                                  ⁢                                          g                      4                                        ⁢                    h                    ⁢                                                                                  ⁢                                          sign                      ⁡                                              (                        g                        )                                                                              -                                      α                    ⁢                                                                                  ⁢                    g                    ⁢                                                                                  ⁢                                          h                      4                                        ⁢                                          sign                      ⁡                                              (                        h                        )                                                                                            )                                                              (                      6            ⁢                          (                              g                -                h                            )                        ⁢                          (                                                g                  2                                -                                  2                  ⁢                                                                          ⁢                  gh                                +                                  h                  2                                            )                                )                                        ⁢                  m        *            =                        (                                    α              ⁢                                                          ⁢                              g                4                            ⁢                                                          ⁢                              sign                ⁡                                  (                  g                  )                                                      -                          α              ⁢                                                          ⁢                              h                4                            ⁢                              sign                ⁡                                  (                  h                  )                                                      -                          2              ⁢                                                          ⁢              α              ⁢                                                          ⁢                              g                3                            ⁢              h              ⁢                                                          ⁢                              sign                ⁡                                  (                  g                  )                                                      +                          2              ⁢                                                          ⁢              α              ⁢                                                          ⁢                              gh                3                            ⁢                              sign                ⁡                                  (                  h                  )                                                              )                          (                                    g              3                        -                          3              ⁢                                                          ⁢                              g                2                            ⁢              h                        +                          3              ⁢                                                          ⁢                              gh                2                                      -                          h              3                                )                    
Finding the Least Squares Model Empirically: Slope-Intercept Model
If the bounded flow range for a pipe segments spans more than a factor of two, then the friction factor may vary significantly over that flow range and there is no analytical expression for the least-squares linear fit of the nonlinear pressure drop relationship. In this case, the preferred approach for developing a least-squares linear fit of the nonlinear pressure drop is a numerical approach.
This approach entails using numerical linear algebra to calculate the value of the slope and intercept using the formula.
      [                            m                                      b                      ]    =                    (                              Q            T                    ⁢          Q                )                    -        1              ⁢          Q      T        ⁢    y  where m is the slope of the line, b is the intercept of the line, Q is a matrix the first column of the matrix Q contains a vector of flow rates ranging from the minimum signed flow rate for the segment to the maximum signed flow rate for the segment, and the second column is a vector of ones.
  Q  =      [                                        q            min                                    1                                      ⋮                          ⋮                                                  q            max                                    1                      ]  
The vector y contains the pressure drop as calculated by the nonlinear pressure drop relationship, at flow rates ranging from the minimum signed flow rate to the maximum signed flow rate. Since the friction factor varies over this flow range, a different value of the nonlinear pressure drop relationship a may be associated with each row of the vector.
  y  =      [                                                      α              min                        ⁢                          q              min                        ⁢                                                        q                min                                                                                        ⋮                                                                α              max                        ⁢                          q              max                        ⁢                                                        q                max                                                                        ]  
As an example, consider the following data from a nonlinear pressure drop model:
Flow,Change in squared pressure,kg/sPa22.07.73.012.14.017.95.025.36.034.17.044.3Given this data,
            q      min        =    2.0    ,            q      max        =    7.0    ,      Q    =          [                                    2.0                                1                                                3.0                                1                                                4.0                                1                                                5.0                                1                                                6.0                                1                                                7.0                                1                              ]        ,            and      ⁢                          ⁢      y        =                  [                                            7.7                                                          12.1                                                          17.9                                                          25.3                                                          34.1                                                          44.3                                      ]            .      Applying the formula
            [                                    m                                                b                              ]        =                            (                                    Q              T                        ⁢            Q                    )                          -          1                    ⁢              Q        T            ⁢      y        ,it is determined that the parameters of the least-squares linear fit are m=7.33 and b=−9.40.
Finding the Least Squares Model Numerically: A Slope Only Model
In some instances, If the flow range includes transition turbulent flow, includes laminar flow, or includes both turbulent and laminar flow regimes, there is no analytical expression for the least-squares linear fit of the nonlinear pressure drop relationship. In this case, the preferred approach for developing a least-squares linear fit of the nonlinear pressure drop is a numerical approach.
This approach involves calculating the value of them=(qTq)−1qty where m is the slope of the line, q is a vector of flow rate values ranging from the minimum signed flow rate for the segment to the maximum signed flow rate for the segment
  q  =      [                                        q            min                                                ⋮                                                  q            max                                ]  
The vector y contains the pressure drop as calculated by the nonlinear pressure drop relationship, at flow rates ranging from the minimum signed flow rate to the maximum signed flow rate. Since the friction factor varies over this flow range, a different value of the nonlinear pressure drop relationship a may be associated with each row of the vector
  y  =            [                                                                  α                min                            ⁢                              q                min                            ⁢                                                                q                  min                                                                                                          ⋮                                                                              α                max                            ⁢                              q                max                            ⁢                                                                q                  max                                                                                        ]        .  
As an example, consider the following data from a nonlinear pressure drop model:
Flow,Change in squared pressure,kg/sPa2−3.0−24.2−2.0−7.5−1.0−1.00.00.01.01.02.07.5Given this data,
            q      min        =    2.0    ,            q      max        =    7.0    ,      q    =          [                                                  -              3.0                                                                          -              2.0                                                                          -              1.0                                                            0.0                                                1.0                                                2.0                              ]        ,            and      ⁢                          ⁢      y        =                  [                                                            -                24.2                                                                                        -                7.5                                                                                        -                1.0                                                                        0.0                                                          1.0                                                          7.5                                      ]            .      Applying the formula m=(qt q)−1qT y, it is determined that the parameter of the least-squares linear fit is m=5.51.
Choosing the Most Appropriate Linear Model
Described herein are several methods for calculating the best linear fit of the nonlinear pressure drop relationship, given the minimum and maximum flow rates in each pipe segment under a range of scenarios including those in which one or more customers is supplied with an increased flow rate. Also, described is how to find the best slope-only linear model, given the minimum and maximum flow rates. An open question is in which situations it is appropriate to use the slope/intercept model, and in which situations it is best to use the slope-only model. A key principle here is that the linear model should always give the correct sign for the pressure drop. In other words, for any linear model exercised over a bounded flow range, the sign of the predicted pressure drop should be consistent with the flow direction. Pressure should decrease in the direction of the flow. Note that the slope-only model has an intercept of zero, and thus the slope-only model will show sign-consistency regardless of the flow range. So, a slope-intercept model should be used unless there is a point in the allowable flow range where there would be a sign inconsistency; if a slope-intercept model would create a sign-inconsistency, then the slope-only model should be used.
Bounding the Error in the Linearized Pressure Predictions for the Pipeline Network
Above described is how to linearize the pressure drop relationship for each pipe in the network by first bounding the range of flow rates which will be encountered in each pipe segment. In accordance with embodiments of the invention, the linearized pressure drop models are used to determine whether it is hydraulically feasible to supply an increased flow of gas to a customer. Although the linearized pressure drop models fit the nonlinear models as well as possible, there will still be some error in the pressure estimates in the network flow solution relative to the pressures that would actually exist in the network given the flows from the network flow solution and the true nonlinear pressure drop relationships. To accommodate this error while still ensuring that pressure constraints are satisfied by the network flow solution, it is necessary to bound the error in the linearized pressure prediction at each node in the network.
To bound the error in the pressure prediction at each node in the network, first, the error in the prediction of the pressure drop for each arc is bound. For pipe arcs, this is done by finding the maximum absolute difference between the linear pressure drop model and the nonlinear pressure drop model in the bounded range of flows for the pipe segment. By definition,
      ps    j    err    =            max                        q          j          min                ≤        q        ≤                  q          j          max                      ⁢                                                              α              j                        ⁢            q            ⁢                                        q                                              -                                    m              j              *                        ⁢            q                    -                      b            j            *                                      ⁢              ∀                  j          ∈                      P            .                              
For control arcs, the maximum error in the prediction of the change in pressure associated with the arc depends on the type of arc. Some control elements, such as valves in parallel with variable speed compressors, have the capability to arbitrarily change the pressure and flow of the fluid within certain ranges, and for these there is no error in the pressure prediction. Other types of control elements, such as nonlinear valves, may be represented by a linear relationship between pressure drop and flow based on the set valve position. For these, there may be a potential linearization error similar to that for pipes. In what follows, it is assumed without loss of generality that psjerr=0 ∀∈C.
Next, a known reference node r in the network is identified. This is a node where the pressure is known with some bounded error. Typically, the reference node is a node which is incident from a pressure control element arc. The maximum absolute pressure error for the reference value can be set to zero, or it can be set to some small value associated with the pressure tracking error associated with the pressure control element.
To compute the error associated with nodes in the network other than the reference node, the undirected graph representing the pipeline network is converted to a weighted graph, where the weight associated with each pipeline arc is the maximum absolute pressure error for the pipe segment. The shortest path is then found, in the weighted graph, between the reference node and any other target node.
In a shortest-path problem, given is a weighted, directed graph G=(N, A), with weight function w: A→R mapping arcs to real-valued weights. The weight of path p=(n0,n1, . . . , nk) is the sum of the weights of its constituent arcs:
      w    ⁡          (      p      )        =            ∑              i        =        1            k        ⁢                  ⁢                  w        ⁡                  (                                    n                              i                -                1                                      ,                          n              i                                )                    .      
The shortest-path weight from n to m is defined by
      δ    ⁡          (              m        ,        n            )        =      {                                                      min              ⁢                                                          ⁢                              {                                                                            w                      ⁡                                              (                        p                        )                                                              ⁢                                          :                                        ⁢                                                                                  ⁢                    m                                    ⁢                                      →                    p                                    ⁢                  n                                }                                                                        if              ⁢                                                          ⁢              there              ⁢                                                          ⁢              is              ⁢                                                          ⁢              a              ⁢                                                          ⁢              path              ⁢                                                          ⁢              from              ⁢                                                          ⁢              m              ⁢                                                          ⁢              to              ⁢                                                          ⁢              n                                                            ∞                                              otherwise              .                                          ,      
A shortest-path from node m to node n is then defined as any path p with weight w(p)=δ(m, n).
In the weighted graph used here, the weight function is the maximum absolute pressure prediction error associated with the pipe segment connecting the two nodes. To compute the shortest-path weight δ(m, n), an implementation of Dijkstra's algorithm can be used (see Ahuja, R. K., Magnanti, T. L., & Orlin, J. B. (1993). Network flows: theory, algorithms, and applications.) The maximum pressure error for the target node is the maximum pressure error for the reference node plus the shortest path distance between the reference node and the target node. In mathematical notation,psmerr=psrerr+δ(r, m)where the weight function for the shortest path is wj=psjerr.
If a pipeline network has more than one pressure reference node r1, . . . ,rn, then one calculates the shortest path between each reference node and every other reference node. The pressure error is then bounded by the minimum of the quantity psrerr+δ(r, m) over all reference nodes:
      ps    m    err    =            min              r        ∈                  {                                    r              1                        ,                                                  ⁢            …            ⁢                                                  ,                          r              n                                }                      ⁢                  {                              ps            r            err                    +                      δ            ⁡                          (                              r                ,                m                            )                                      }            .      
Determining Whether It is Hydraulically Feasible to Supply an Increased Flow Rate to a Customer
Described herein is a method for analyzing a scenario where a customer receives additional product by 1) bounding the minimum and maximum flow rate for each pipe segment in a computationally efficient fashion; 2) computing an accurate linear approximation of the nonlinear pressure drop relationship given the bounded flow range; 3) bounding the pressure prediction error associated with the linear approximation. Now it can be determined whether it is hydraulically feasible to supply additional product to a customer, that is, to determine whether there is an increased flow scenario which 1) satisfies constraints associated with the conservation of mass and momentum; 2) is consistent with bounds on the flow delivered to each customer, 3) satisfies pipeline pressure constraints with appropriate margin to accommodate errors associated with the linearization of the nonlinear pressure drop relationship. The governing equations are summarized here. Node mass balance
The node mass balance stipulates that the total mass flow leaving a particular node is equal to the total mass flow entering that node.
            d      n        +                  ∑                  j          ❘                                    (                              n                ,                j                            )                        ∈                          A              in                                          ⁢                          ⁢              q        j              =                    ∑                  j          ❘                                    (                              n                ,                j                            )                        ∈                          A              out                                          ⁢              q        j              +          s      n      
Node Pressure Continuity
The node pressure continuity equations require that the pressure of all pipes connected to a node should be the same as the pressure of the node.psjin=psnnode∀(n, j)∈Ain psjout=psnnode∀(n, j)∈Aout 
Linearized Pressure Drop Mode
We have shown how to develop a linear pressure drop model of the formpsjin−psnout=mjqj+bj.
Pressure Constraints at Nodes
At nodes in the pipeline network, there are minimum and maximum pressure constraints. These constraints must be satisfied with sufficient margin, namely psnerr, to allow for potential inaccuracy associated with the linearized pressure drop relationships:psnmin+psnerr≤psnnode≤psnmax−psnerr, ∀n ∈N . 
This ensures that the pressures constraints will be satisfied even when the nonlinear pressure drop model is used to calculate network pressures based on the flow values associated with the network flow solution. Above, we have shown how to compute psnerr using Dijkstra's algorithm for a certain weighted graph.
Production Constraints
This constraint specifies the minimum and maximum production rate for each of the plants.snmin<sn<snmax
The governing equations can be combined to formulate the following linear program to determine the maximum flow rate of a gas that can be supplied to customer k.
GIVENdncurr ∀n ∈ NCurrent customer demand rate in node n(mj, bj) ∀ j ∈ PLinearized pressure drop model for pipe jpsnerr ∀ n ∈ NMaximum squared pressure error fornode n, given linearized pressure dropmodelssnmin < sn < snmaxMinimum and maximum production ratesat node nCALCULATEqj ∀ j ∈ AFlow rate in arcssn ∀ n ∈ SProduction rate in supply nodedn∀ n ∈ DUpdated rate supplied to customerspsnnode ∀ n ∈ NSquared pressure at each nodepsje ∀ j ∈ ASquared pressure at the ends of each arcIN ORDER TO MAXIMIZEdkMaximum flow rate of gas which can besupplied to customer of interestSUCH THATdn + Σj|(n, j)∈Ain qj = Σj|(n, j)∈Aout qj + sn ∀ n ∈ NNode mass balancepsjin = psnnode ∀(n, j) ∈ AinNode pressure equality constraintspsjout = psnnode ∀(n, j) ∈ AoutNode pressure equality constraintspsjin − psjout = mj qj + bj ∀ j ∈ PLinearized pressure drop model for pipespsnmin + psnerr ≤ psnnode ≤Pressure bounds with margin for errorpsnmax − psnerr, ∀n ∈ Nsnmin < sn< snmax ∀n ∈ SProduction boundsdncurr < dn < dncurr + Δn ∀n ∈ DDemand bounds when one morecustomers accept additional product
The above linear program can be quickly solved by a wide variety of linear programming solvers, including those in MATLAB, Gurobi, and CPLEX. Note that additional linear constraints, such as min or max flow rates in certain arcs, can easily be added to the above linear program. The primary result of the linear program, dk , is the maximum flow rate to customer k that is hydraulically feasible. If this amount is significantly greater than the current flow rate of gas being supplied to the customer, then it is hydraulically feasible to offer an increased flow of gas to the customer.
Estimating Whether a Customer Has Latent Demand for A Gas
Above described is a computationally efficient method to determine whether it is hydraulically feasible to supply an increased flow rate of a gas to a customer. In order to increase the capacity factor of a gas pipeline network, it is also important to determine whether a customer has latent demand for the gas. In order to efficiently increase the capacity factor of the gas pipeline network, it is important to have a means to automatically determine whether latent demand exists on a frequent and regular interval, without the need to query the customer.
Embodiments of the invention incorporate a classification tree which uses intrinsic and extrinsic factors to determine whether latent demand for the gas exists. A classification tree is a machine learning construct which uses certain features to predict an outcome. The classification tree can be represented as a binary tree, where the classification tree starts at the root of the tree. A series of binary decisions are made, based on the values of intrinsic and extrinsic factors available to the operator of the gas pipeline network.
A classification tree for whether there is latent demand for a gas may be generated automatically using historical data as to whether the customer accepted an increased flow rate of gas when it was offered to them.
Other machine learning models may be used within the scope of the present invention, such as a support vector machine.
Consider the example where the gas pipeline network is for the production, transmission, and distribution of hydrogen, and a customer for hydrogen gas is a petroleum refinery. The operator of the pipeline network may have information on intrinsic and extrinsic factors associated with latent demand for hydrogen. Examples of intrinsic factors affecting the latent demand for hydrogen are: 1) the change in consumption of hydrogen by the refinery over the past three hours, 2) time-of-day, 3) the crude slate for a particular refinery, 4) day of week, and 5) time since last call. Examples of extrinsic factors affecting the latent demand for hydrogen include 1) the retail price of gasoline, 2) the price of ultra-low diesel sulfur, 3) the rate of imports of petroleum for the region in which the refinery is located, 4) the price of natural gas, and 5) the spread between the prices of sweet and sour crude.
FIG. 8 is an example of a classification tree that might be used to determine whether or not a petroleum refinery has latent demand for hydrogen.
A variety of machine learning techniques, other than classification trees, may be used to estimate whether a customer has latent demand for an industrial gas. Other techniques that might be used include logistic regression, linear discriminant analysis, Fisher discriminant analysis, and support vector machines.
Receiving an Updated Request Rate for the Gas
If it is hydraulically feasible to supply an increased flow rate of gas to the customer, and it is estimated that the customer has latent demand for the gas, then an updated request rate is received from the customer. The updated request rate may be received telephonically, by email, or by other electronic means. Often, the updated request rate would be provided in response to an offer from the operator of the industrial gas network to provide an increased flow of the gas. Typically, the updated request rate would be for a certain flow rate of the gas which is as much as Δn units greater than the current rate. In describing how a network flow solution is calculated below, the newly updated request rate for customer n is represented by the variable dn.
Calculating a Network Flow Solution Using an Updated Customer Request
After it has been determined that it is hydraulically feasible to provide an increased flow of the gas to a customer, and it has been estimated that the customer may have latent demand for the gas, and an updated request rate has been received, embodiments of the invention calculate a new network slow solution. The network flow solution is calculated using the linearized pressure drop models that were described above. Embodiments of the invention use a linear program as follows:
GIVENdn∀ n ∈ DUpdated rate to be supplied to customers(mj, bj) ∀ j ∈ PLinearized pressure drop model for pipe jpsnerr ∀ n ∈ NMaximum squared pressure error fornode n, given linearized pressure dropmodelssnmin < sn < snmaxMinimum and maximum production ratesat node nCALCULATEqj ∀ j ∈ AFlow rate in arcssn ∀ n ∈ SProduction rate in supply nodepsnnode ∀ n ∈ NSquared pressure at each nodepsje ∀ j ∈ ASquared pressure at the ends of each arcSUCH THATdn + Σj|(n, j)∈Ain qj = Σj|(n, j)∈Aout qj + sn ∀∈ NNode mass balancepsjin = psnnode ∀(n, j) ∈ AinNode pressure equality constraintspsjout = psnnode ∀(n, j) ∈ AoutNode pressure equality constraintspsjin − psjout = mj qj + bj ∀ j ∈ PLinearized pressure drop model for pipespsnmin + psnerr ≤ psnnode ≤ psnmax − psnerr, ∀n ∈ NPressure bounds with margin for error
The linear program may be solved using any of a variety of linear program solvers, including those found in Matlab, CPLEX, or Gurobi.
Controlling the Gas Pipeline Network Using the Network Flow Solution
Once the network flow solution has been computed, it can be used to control the gas pipeline network. Flow control elements receive setpoints which are identified using the network flow solution.
There are two representations of control elements in the undirected graph representation of the network. First, nodes associated with supply or demand are control elements, and the network flow solution indicates the supply or demand flow that should be associated with each plant or customer in the network. Second, in some networks there are also control arcs (representing compressors, valves, or a combination of compressors in valves). The network flow solution indicates the flows and pressures that should be accomplished by these control elements.